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Abstract 

In this paper we discuss how partial knowledge of the density of states 
for a model can be used to give good approximations of the energy dis- 
tributions in a given temperature range. From these distributions one 
can then obtain the statistical moments corresponding to eg the internal 
energy and the specific heat. These questions have gained interest apro- 
pos of several recent methods for estimating the density of states of spin 
models. 

As a worked example we finally apply these methods to the 3-state 
Potts model for cubic lattices of linear order up to 128. We give estimates 
of eg latent heat and critical temperature, as well as the microcanonical 
properties of interest. 

1 Introduction 

When studying a statistical mechanical model the most complete information 
is given by the density of states function. From complete knowledge of the 
density of states one can immediately work with the microcanonical ensemble 
and of course also compute the partition function and through it have access 
to the canonical ensemble as well. The main problem here is that computing 
the density of states for systems of even very modest size is typically very hard. 
However, recently several sampling schemes which strive to approximate the 
density of states have appeared. One recent metho d was given [WLOlj and in 
[WS02| several such methods were given, and in H RA + 04b| all of the later 



methods as well as several others were united in a common framework. 

For work in the microcanonical ensemble the mentioned methods give all 
the information needed. Using them one can find the density of states in an 
energy interval around the critical region and that is all that is needed for 
most investigations of the critical properties of the model. The microcanonical 
ensemble is more refined than the canonical ensemble in that every equilibrium 
measure for the canonical ensemble is found among the equilibrium measures 
for the microcanonical ensemble, but for some models there are microcanonical 
equilibria which are not present in the canonical ensemble. For a fuller survey of 
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the mathematical theory of ensemble equivalence see |TET04j and its references. 
This means that all properties of the thermodynamic limit can be obtained via 
the microcanonical ensemble. 

However, even in view of what has been said the canonical ensemble has its 
own interest for finite systems. Among other things it governs the behaviour of 
many sampling algorithms and for systems where we have ensemble nonequiva- 
lence its dynamic can be very interesting. In order to reconstruct the canonical 
ensemble one would in principle need to know the density of states for all values 
of the energy E. However, using methods as in |HRA+04bj this is very costly, 
and also not needed for work in the microcanonical ensemble. 

Our aim is to look at how density of states data from a restricted interval 
of energies can be used to get an approximation of the energy distribution of 
the canonical ensemble for some range of couplings K. Thanks to the strong 
concentration of the energy distributions we will see that one can obtain a very 
good approximation of the energy distribution and through its moments most of 
the standard thermodynamical properties. This will be demonstrated first in a 
case where we know the exact partition function, the Ising model on the 256 x 256 
square lattice, and then for a case where we have ensemble nonequi valence: the 
3-state potts model on the 3-dimensional cubic lattice. All in all we find that 
with data collected with the methods of [HRA^t04b] in mind one can get a good 
picture of the canonical ensemble as well as the microcanonical. In fact, thanks 
to knowing the density of states for a full interval of energies we will be able 
to reconstruct the canonical ensemble for all couplings in some interval rather 
than just those used in the sampling process. 



2 Notation 

Let us define what we need in terms of the Ising model. Later on, when the 
Potts model is our subject, we will redefine some quantities, but our general 
discussion will be held in terms of the Ising model. Let G be a graph on n 
vertices V = {1, . . . , n} and m edges. A state is a function s : V — > Q where 
Q = { + 1,-1} and we say that vertex i has spin Sj. The energy of a state is 
defined as E(s) = Y]^ SiSj where the sum is taken over all edges ij of the graph 
and we have —m < E < m. The magnetisation of s is defined as M(s) = Sj 
so that — n < M < n. 

A normalised energy and magnetisation will often be used, here defined as 
U = E/m and /j, — M/n so that —1 < U, n < 1. The number of states 
having energy E and magnetisation M is denoted a(E,M). The number of 
states at energy E, or, the density of states, is denoted a(E), where, of course, 
a(E)=^2 M a(E,M). 

From quotients of a(E) we obtain what we will refer to as the coupling 
function 

where U = E/m and £ is the difference between two consecutive energies. The 
very fundamental entropy function 

S(U) = 



2 



is of course related to the coupling function through 

77 

K{U) = S'(U) 

m 

See [HRA + 04b] for proofs and further details. 

The partition function is defined for all graphs as 

Z(K,H) = a(E,M) exp(K E + H M) 

E.M 

where K and H are the dimensionless coupling and external field respectively. 
When the external field is zero we simplify as 

Z{K) =J2a(E) exp(KE) 

E 

As a convention we will write our coupling dependent quantities in a calligraphic 
font, such as Z(K). 

The central moments of a random variable X are defined 

Gi = ((x- (x)Y), i = o,i,... 

where <7q = 1, a% = and <72 = Var(X). The standard deviation is written 
a = ^Jo2- The cumulants Ki of a distribution, or, the i:th derivatives of log-E, 
can be expressed in terms of moments. For the first few we have 



--Ki = (E) 

-~k 2 = Var (E) = o- 2 



=«3 = ^3 



o 2 
-Ki = <J& — 6 
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d log Z (K) 

dK 
d 2 log Z(K) 

dK 2 
d 3 log Z{K) 

dK 3 
<9 4 log Z{K) 
dK 1 

The free energy is here defined as 

T(K) - - log Z(K) 
n 

and the reader should note that we have used a simplified version compared to 
its traditional form. The internal energy, specific heat and coupling dependent 
entropy are given by 

m oK m 

C(J0 _ I 

m oK z m 

777 

S{K)=T{K) KU{K) 

n 

We would also like to study the higher derivatives in the form of skewness 

d 3 log Z(K)/dK 3 cr 3 k 3 



(d 2 logZ(K)/dK 2 ) 3/2 (a 2 ) 3/2 ( K2 ) 3/2 
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and (excess) kurtosis 



d 4 log Z(K)/dK 4 _ <74 _ 3 _ K4 
(d 2 log Z(K)/dK 2 ) 2 ~4 



Note that for normal distributions they both evaluate to zero. 

Derivatives with respect to the field H are of course obtained analogously. 
The magnetisation and susceptibility are defined respectively as 

However, what one usually want is the spontaneous magnetisation and suscep- 
tibility. As finite size approximations of these we use 

fi(K) = -{\M\) 
n 

x(J0 = ivar(|M|) 

and assume that these converge to the appropriate limits. 

Given a lattice of side L with L 3 vertices we call L the linear order of he 
lattice. When necessary we will subscript the functions with the linear, as in 
Z L . 



3 Distributions of energy 

In this section we will look at how the distribution of energies for a given coupling 
Kq can be reconstructed. The process is rather straightforward and follows more 
or less by definition, but we will derive it in some detail. We will first derive an 
exact expression for Pr(£ l ) when S'(U) is exactly known, next we look at what 
can be done when only partial knowledge of K(U) is available, and finally we 
consider precision issues for such incomplete reconstructions. 



3.1 From coupling to distribution 

Our first aim is to express Pr(E) in terms of the values of K{U) in some interval 
of energies u < U < v. 

We assume the Boltzmann distribution for the states, that is, if we sample 
at a coupling K the probability for our system being in state s is 

cxp{K E{s)) 
(J ~ Z(K ) 

and consequently the probability for our system being in a state of energy E is 

Pr(E) .~^ 

Recall that we defined a(E) — exp(nS(U)). Then we obtain 

cxp(nS(U)+mUK ) 

Pr(E) = m^) (2) 
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By definition we also have 



pU nU pU rU 

S{U)= / S'(x)dx= / S'{x)dx+ / S\x)dx = A u / K(x) dx 

J -l J-l Ju n Ju 

and trivially 



U = u+ I Idx 



Plugging these identities into Equation [5] and applying only a modicum of alge- 
braic manipulation it simplifies finally into 



Pr (E) = C u exp J K - K{x) dx 



(3) 



Since the outcome is a probability function the constant C u can also be defined 
by normalising so that 

£Pr(£) = l (4) 

E 

We will consider C u further in the next section. Finally, we note in passing 
that the derivative of the probability function with respect to U is m (Kq — 
K(U)) Pr (E). Thus the points where the sign of the derivative changes is de- 
termined by when Kq = K(U). 

Note that we have only defined the function K(U) at discrete points U = 
E/m so we should be somewhat careful with how the integral is taken. If a 
function f{x) is defined at a = xq < x\ < . . . < x p — b then we use a left-point 
rule for integration 



pb p- 1 

/ f(x) dx = f( x i) (^+i 

Ja i=0 



Having reconstructed the distribution of energies the moments and cumu- 
lants are easily retrieved. First the average 

(E) =Y,E?r{E) 

E 

and then the central moments 

E 

and from these we obtain the sought-after estimates of the derivatives by eval- 
uating the cumulants of the distribution so that, for example 



Hi = (74 — 3 a 



2 _ a 4 iogz 

2 



Let us address the issue of derivatives with respect to the field as well. If we 
during our sampling process remembered to collect data on the magnetisation as 
well, then we can reconstruct the spontaneous magnetisation and susceptibility 



5 



as well. Our program should then collect raw moments on the form (\M\ l | E). 
Then the following holds 

fi(K ) = \ {\M\) = \ «|M| \E)) = \ "£{\M\ | E) Pr(E) (5) 

E 

x(Xo) = ivar(|M|) = i ^(|M| 2 | E) Pr{E) - |£>Pr(£^ 

(6) 

The following is a nice alternative way of writing the variance 

Var(|M|) = (Var(|M| | E)) + Var «|Af| | -E 1 )) 

that is, the variance is the sum of the expectation of the variances and the 
variance of the expectations. 



3.2 Reconstruction from incomplete data. 

In the previous subsection we assumed that K(U) was exactly known for all en- 
ergies. When this is the case we have seen that Pr(E) can be exactly determined, 
as one would expect. Let us now assume that our data contains information 
on quotients of consecutive density of states, or rather, that we have available 
estimates for the coupling function K(U) for an interval of energies u < U < v. 

In this situation we can no longer determine Pr(iS) completely from Equa- 
tion [3] since we can no longer compute C„, and we may have errors in our 
estimate for K{E). However, if the expected energy lies well within the interval 
u < U < v we can hope that a good enough approximation can be found, due 
to the rapid decay of the tails of the energy distribution. 

The simplest way to find such an approximation is simply to let Equation 0] 
define an approximate value C' u via 

u<E /m<v 

where Pr u (E) is obtained by using C' u instead of C u in Equation [31 Equation 
[7] is simply a linear equation for C' u and so can be easily solved. Once C' u has 
been found we can use C' u and our estimated K{U) in Equation [3] to compute 
an approximate distribution function for the canonical ensemble. 



3.3 Accuracy of the reconstructed distributions 

Since the canonical ensemble is always determined by the density of states we 
only have two sources of errors: the precision of the original data and the 
truncation error due to not having data from all energies. 

In a perfect world the collected data comes from the entire interval of energies 
— 1 < U < 1. However, normally it suffices for the interval to be wide enough 
to cover the energies at coupling Kq with a high probability. In short, the 
distribution of energies corresponding to Kq must stay in the interval [u, v] with 
a probability close to 1. If [it, v] only covers say, 99% or less of the energies you 
see at Kq, the normalising step in Equation U will produce erroneous results. 
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Given a coupling Kq that is close to the critical coupling K c we expect the 
distribution to be anything but normal (ie gaussian). But, as we move away from 
K c the distribution typically becomes close to normal. For example, at K = 
the distribution is clearly approaching a normal one with increasing system size. 
It has been shown, see |ML73j , that this also holds for Ising systems when K is 
greater than some K\ > K c . Since our approach is somewhat pragmatic we will 
only assume that far enough from K c the energy distributions can be treated 
as roughly normal. 

For how large or small values of Kq does the procedure return a credible 
distribution on [u,v]l Treating the distribution as roughly normal it should be 
enough, for all practical purposes, to make sure that the end-points are at least 
4, if possible 5 and preferably 6 standard deviations a away. 

The probability density of a normally distributed variable is 



fix) 



/2vr 



cxp(-a; 2 /2) 



We will translate it slightly to the right, ie take f(x — p) for p > 0, and cut it 
off at x = 0. Now let X be a random variable with probability density 



x < o 







A(p) 



where A(p) is the mass of probability on x > 0, ie 

/•OO 

A(p)= / f(x~p)dx 



so that g(x) becomes a cut-off, but otherwise normal looking, probability density 
on the real axis. For which p is the y-axis, ie the cut-off point, located k standard 
deviations a away? Numerical calculations gave Table [T] below and in Figure [T] 
the distribution functions are shown for k = 2, 3, 4. In the Table we also list the 
errors 

£i = Wi(g) - Ki(/)| 

that is, we take the difference in cumulants for the cut-off density g and the 
normal density / translated k standard deviations. These errors are of course 
very idealised, being based on normal distributions, and should be considered 
rough guidelines. For any particular distribution we will see different errors and 
especially the higher cumulants will deviate from these. 
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Table 1 : Peak location p and cumulant errors for a cut-off normal distribution 
with (X) = ka. 
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Figure 1: Cut-off normal distributions such that (X) = k a for k = 2, 3, 4. 

4 The 2-dimensional Ising model 

We will employ the 2D Ising model as a test bed for our method. Recall that 
the critical coupling is K c — arct anh (\/2 — 1) w 0.4407 and that the critical 
energy is U c = \ w 0.7071. In |HRA + 04a] we computed the exact partition 
function for the 256 x 256-lattice with periodic boundary. However, since the 
largest density of states a(0) has 19726 digits we will take the liberty of doing 
all actual computations with 50 digits numerical precision instead. 

Suppose now that we have collected data on K{U) for u = 0.6 < U < 0.8 = 
v, an interval comprising 6554 energies. In Figure[2]we plot K{U) and K'(U) for 
L = 256. From the exact (50 digits) coupling function on the interval [0.6,0.8] 




6 0.65 0.7 0.75 0.8 0.6 0.65 0.7 0.75 0.8 



Figure 2: K(U) and K'(U) for 256 x 256-lattice. 

we reconstruct the distribution, ie the probability density function, of energies 
at K c using Equation [3] Let denote the relative error of the z:th cumulant 
where we compare the cumulant Ki of the reconstructed energy distribution with 
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the i:th derivative of logiJ, ie 



d l log Z/dK l 

The relative errors ei, £ 2 , £3, £4 are negligibly small, less than 1-10 -30 . However, 
this distribution lives clearly in the middle of our energy interval [0.6,0.8], the 
lower bound being 14cr below and the upper bound 12a above the mean. In 
Table [2] we compute relative errors of the cumulants when the coupling cor- 
responds to a cut-off distribution with (E) located k a from the lower bound 
u = 0.6 for k = 2, 3, 4, 5, 6 and in Table [3] we do the corresponding at the other 
end of the interval so that (E) is ka from the upper bound v = 0.8. Figure [3] 
shows the probability densities at K = 0.423780, K = K c and K = 0.454942. 
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Table 2: Relative errors of cumulants for cut-off distribution with (X) = u + ka. 
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Table 3: Relative errors of cumulants for cut-off distribution with (X) = v — ka. 



We also computed the cumulant errors at K c with the upper and lower bound 
of the energy interval located 6 a away. For L — 32, 64, 128, 256 the errors are 
quite small, e x < MO" 12 , e 2 < 2-10" 10 , e 3 < 2-10~ 8 and e 4 < 6-10~ 8 . For 
L < 16, the errors become larger, but on the other hand for such small graphs 
it is easy to collect a complete set of K (C/)-data instead of only a short interval. 

Regarding the magnetisation and susceptibility we have no way of comparing 
the reconstructed values with exact values. We have simply run the Metropolis 
method at 10 different temperatures in the vicinity of K c (0.42 < K < 0.46) and 
collected magnetisation moments at each energy level. Using these data and the 
exact if-function we can reconstruct the magnetisation and susceptibility at any 
temperature in that region using Equation [5] and [5] adding data as prescribed 
in HRA + 04~b] , See Figure [3] for a snap-shot of the result. The reconstructed 



curves agrees very well with the sampled data. 
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Figure 3: Probability densities of energies at K\ = 0.423780 (left), K c (middle) 
and K 2 = 0.454942 (right) for 256 x 256-lattice. At K x and K 2 the interval 
bounds are 5a away. Probability Pr (E) on the y-axis and energy U = E/m on 
the x-axis. 
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Figure 4: fi(K) and x(K) for 256 x 256-lattice. 



5 The free energy 

By definition we have that 

m f K 

T{K) = F(0) + — J U{x)dx 

where the constant F(0) = log 2 for the Ising model, and .F(O) = logg for the 
g-state Potts model. Having evaluated U (K) for a number of values of K this 
is of course easily accomplished. Unfortunately, this formulation implies that 
we have collected data so that the energy distribution can be reconstructed for 
K > 0. For smaller graphs this is of course perfectly alright but for large graphs 
this was exactly what we wanted to avoid. However, due to the well-behaved 
nature of the internal energy U(K) we can circumvent this problem. Suppose 
we have reconstructed the internal energy for two system sizes L\ and L 2 , where 
L\ < L 2 - Suppose further that we have Ul 1 (K) for < K < b\ and Ul 2 IK) for 
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a < K < 6 2 where < a < K c < 62 < &i- Then, for a < if < 62 we have 

ml" 1 m T K 

T(K) = + - / W Ll (x) dx + - U L2 (x) dx + e 

"JO " J a 

where e is an error term. 

How big is the error? Let / and g be continuous functions on the interval 
[a, c] with a < b < c. Then the following elementary calculation gives an 
estimate: 

f-C r-b f-C 

I g(x) dx — g(x) dx + g(x) dx = 

J a J a Jb 

= / g(x) - f(x) + f(x)dx+ / g{x)dx = 

J a J b 

pb pc 

= / f(x) dx + g(x) dx + e 

Ja Jb 

where e is the error term 

e = / g(x) - f(x)dx 

J a 

which gives the very simple but useful estimate 

|e|<(6-a) max \g{x)-f(x)\ (8) 

a<x<o 

Since the internal energy function is an increasing and, in fact, convex function, 
it is easy to establish the maximum. The integration is numerical so it is impor- 
tant to evaluate IA (K) at points chosen densely enough, with special attention 
to values close to K c where U (K) is expected to change rapidly. 

5.1 A worked example for the 2D Ising model 

Here our goal is to compute the free energy at K = K c for the 256 x 256 2D Ising 
model by using a sequence of system sizes, L = 32, 64, 128, 256, and formulate 
the method as 

0.15 /-0.30 

+ 



f-V.10 /-U.dU 

^256 {K c ) = log 2 + 2 / U 32 (x) dx + 2 U 6i [x) dx 

JO JO. 15 

/•0.40 rK a 

+2 / U 128 (x)dx + 2 / W 2 56 (x) dx 

JO. 30 JO. 40 



though the integration will of course be numerical. To evaluate the Ul{x) at 
the couplings indicated by the integral boundaries we use the exact (to 50 digits 
precision) if -functions at intervals wide enough to keep the endpoints 6 a away 
for each energy distribution. Picking simple values gives the following energy 
intervals, coupling intervals and step lengths used for evaluating Ul(K). 

L = 32, -0.15 < U < 0.30, 0.00 < K < 0.15, h = 0.0010 

L = 64, 0.05 < U < 0.45, 0.15 < K < 0.30, h = 0.0005 

L = 128, 0.30 < U < 0.60, 0.30 < K < 0.40, h = 0.0002 

L = 256, 0.50 < U < 0.75, 0.40 < K < 0.441, h = 0.0001 
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The intervals for K were chosen to give good overlaps, and are of course model 
dependent. 

Numerical integration of these data points, using eg the trapezoidal method 
gives that T 256 (K C ) « 0.92970521. The correct answer is T 2 ^{K C ) = 0.92970516 . 
giving an error of 5.5 T0~ 8 . 

The error contribution given by Equation [5] is only of the order of 6T0~ 12 . 
The main error source is actually the numerical integration. As is well-known, 
numerical evaluation of J f(x) dx with the trapezoidal rule gives the error term 

(b-a)h 2 
e = - 12 / (0. a<ti<b 

where h is the step length. Since the function we integrate is U(K) its second 
derivative is 

1 dHogZ(K) k 3 

U \ K > = = _ 

and it is at a very little extra cost we evaluate the third cumulant when we 
already have the distribution. 

For example, the error contributed from 

6=0.15 

U 32 (x) dx 



a=0 



is at most 



2(b-a)h 2 , . 2 0.15- 0.001 2 „„ A „ H 
— : - — max Uon(x) = — y - '- 1.67 « 4.2-10" 8 

12 a<x<b 32K ' 12 

and the errors contributed by the other integrals are at most respectively 3.410 -8 
for L = 64, 1.6- 10~ 8 for L = 128 and 5.8- 10~ 8 for L = 256 and they sum up to 
1.5 TO" 7 which is clearly larger than the actual error we received. 



6 An example with a first order phase transi- 
tion: The 3-dimensional 3-state Potts model 

For this model we need to redefine some of our quantities. A state is here a 
function s : V — * Q where Q is a set of q distinct elements, eg Q — {1, . . . , q}. 
The energy is defined as E(s) = <5(s£, sj), where 5{x,y) is the Kronecker- 
delta, so that < E < m. We normalise as before and let U = E/m so that 
< U < 1. Let rjj — ^ S(si,j), ie the number of vertices having spin j. For the 
Potts model, the definition of magnetisation M varies slightly in the literature. 
For example, M = max(?7i, . . . , rj q ) is sometimes used, but M = r/f + . . . + r/ 2 
is the one used here. This means that n 2 /q < M < n 2 and we normalise by 
taking jl — /i = M/n 2 so that 1/q < p, < 1. Having defined these quantities 
their physical versions follow accordingly. 



6.1 The sampled data 

The data w ere generate d and collected by using the sampling method described 
in detail in H RA + 04b and we refer to that paper for further details. Since this 
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model is conjectured to have a first order phase transition and cluster methods 
thus are expected to have exponential mixing time |CJF + 99 we opted for a 
highly optimised single spin Metropolis method. We used up to a few hundred 
independent spin systems, which after slowly being brought to the right coupling 
were given a few days or weeks, depending on their size, of continuous running 
for mixing. The length of the sampling runs were of the same order. The 
sampled also passed the consistency test from [HRA + 04b] . 

For the smaller lattices we collected at least 10000 measurements per energy 
level, often orders of magnitude more. For the larger lattices, say L > 32, 
this quickly becomes difficult. For L = 96, 128 we did not manage to fill out all 
energy levels inside the energy jump at critical the coupling, though these empty 
levels are very few relatively speaking. From these data we then constructed 
the coupling function K{U). 

The coupling functions K(U), from which the energy distributions are gen- 
erated, are shown in the left plot of Figure and the right plot shows the 
magnetisations n(U). Note that the if- functions behave rather different from 
that of the Ising model in Figure [2J Here the if-functions have their own set of 
critical points and they are listed in Table HJ Let U~ and U + be the locations 
of the maximum and the minimum respectively of K(U). The corresponding 
values of K at these points are denoted K~ and K + respectively. We define the 
latent heat as = U + — U~ . Let also /i + and \i~ denote the magnetisation 
at respectively U + and U~ . 




0.53 0.54 0.55 0.56 0.57 0.58 0.59 



Figure 5: K(U) (left) and n(U) (right) for L > 16. 



L 


U~ 


U+ 


U± 


K~ 


K+ 


i l 


M + 


6 


0.54502 


0.60677 


0.06174 


0.553398 


0.548924 


0.4174 


0.5125 


8 


0.54138 


0.59245 


0.05101 


0.553159 


0.549088 


0.3930 


0.4779 


12 


0.53851 


0.57979 


0.04128 


0.552756 


0.549367 


0.3694 


0.4442 


16 


0.53721 


0.57327 


0.03606 


0.552507 


0.549588 


0.3580 


0.4264 


24 


0.53603 


0.56884 


0.03281 


0.552133 


0.549861 


0.3474 


0.4119 


32 


0.53564 


0.56717 


0.03153 


0.551905 


0.550012 


0.3427 


0.4057 


48 


0.53478 


0.56645 


0.03167 


0.551625 


0.550151 


0.3380 


0.4012 


64 


0.53457 


0.56796 


0.03340 


0.551465 


0.550233 


0.3363 


0.4029 


96 


0.53337 


0.56816 


0.03479 


0.551242 


0.550284 


0.3344 


0.4019 


128 


0.53090 


0.56785 


0.03694 


0.550826 


0.550306 


0.3336 


0.4004 



Table 4: Critical points and values of the combinatorial functions. 
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6.2 The reconstructed ensemble 



Next we used the {/-dependent functions from the previous section to reconstruct 
the if-dependent quantities in the way described earlier in the paper. Let us 
here stress that none of the functions plotted here were sampled directly, ie we 
did not keep track of the variance and expectation of the energy in the sampling 
runs and everything here is based on the microcanoncial data. 

First we will show a quick gallery of pictures of the physical, ie coupling 
dependent, quantities that were defined in Section [2] In Figure El we show the 
free energy T and the internal energy 11 in a narrow region around K c for the 
larger lattices. The dramatic jump in the energy of course leads to a similar 
behaviour in the entropy S(K) shown in Figure [7] This is also seen in the 
magnetisation fi(K) in the same figure. We find that everything agrees well 
with the expected first order nature of the phase transition. 




.5503 0.5504 0.5505 0.5506 0.5507 



0.5503 0.5504 0.5505 0.5506 0.5507 



Figure 6: Free energy T(K) (left) and internal energy U(K) (right) for L > 32. 




0.5503 0.5504 0.5505 0.55 



Figure 7: Entropy S(K) (left) and magnetisation fi(K) (right) for L > 32. 



The specific heat C(K) is shown for L = 64 in the left plot of Figure [HI The 
maximum of this quantity grows very fast with L and to be able to compare 
them for several L the right plot shows their logarithm. This can also be said 
about the skewness and kurtosis in Figure [U These quantities changes sign in a 
number of critical points so taking logarithms is not advisable. The distributions 
go through a sharply bimodal phase as the coupling moves past K c . Define K* 
as the point where the specific heat has its maximum. What do the distributions 
at this point look like? In the left plot of Figure [10] the probability densities 
p(x) of the normalised variable x = (E — (E))/a(E) at K* are shown, while the 
right plot shows the distribution functions (accumulated densities) defined as 



<S>(x) 



p{t) dt. 
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0.5503 0.5504 0.5505 0.5506 0.5507 



Figure 8: Specific heat C(K) for L = 64 (left) and its logarithm for L > 32 
(right). 




Figure 9: Skewness Ti(K) (left) and kurtosis T 2 (K) (right) for L = 64. 



Note, by the way, that the peaks in the probability densities are very close ±ct. 
In Tabled] the data connected with K* are listed. 




-2 -1.5 -1 -0.5 0.5 1 1.5 2-2 -1.5 -1 -0.5 0.5 1 1.5 2 



Figure 10: Normalised probability densities (left) and distribution functions 
(right) vs cr at K* for L > 32. 



6.3 Asymptotics 

In this section we will see how some of the values in the tables above scale 
with the linear order. First we wish to establish the critical coupling K c . We 
have three separate sequences of critical points which should all converge to 
K c , namely K* , K + and K~ . The coefficients of the fits described below are 
collected in Table 16.31 

The sequences K + and K~ from Tabled] have the nice feature that they are 
monotone; K + is increasing and K~ is decreasing. Unfortunately though, the 



15 



L 




K* 


C(K*) 


U(K*) 




S(K*) 




6 


o 


555045 

'.J "7 '.J V 7 1 "7 


5.1279 


63766 


I 


78444 


72265 


0.5557 


8 


o 


552821 


7.2889 


61575 

\ > ■ V7 1 • / J 7 


I 


77697 


0.75577 


0.5130 


12 


o 


551143 


12.539 


59190 


1 


77220 


79355 


0.4643 


16 





550665 


19.584 


0.57983 


1 


77093 


0.81306 


0.4381 


24 





550460 


41.929 


0.56831 


1 


77034 


0.83185 


0.4118 


32 





550462 


81.880 


0.56318 


1 


77028 


0.84026 


0.3996 


48 





550502 


255.97 


0.55921 


1 


77032 


0.84678 


0.3901 


64 





550530 


614.58 


0.55821 


1 


77036 


0.84842 


0.3877 


96 





550432 


1790.0 


0.55550 


1 


77020 


0.85291 


0.3821 


128 





550353 


3422.7 


0.55279 


1 


77008 


0.85738 


0.3763 



Table 5: Critical points K* and values of the physical quantities. 



i4f - -sequence appears slightly blemished for L = 128, as is the [/"-sequence. 
Even so, after discarding that particular point and assuming that the sequences 
stay monotone also for larger L we then have upper and lower bounds on K c . 
Then 0.550306 < K c < 0.551242, a rather wide interval. The story is different 
for the if* -sequence, it sometimes increases, sometimes decreases, but we see 
nothing amiss with the value for L — 128. 

To establish a K c we have attempted a simple fit of the form Co + c% L~ x for 
some coefficients Co, c\ and exponent A to the data for K*. To find the param- 
eters we used Mathematica's non-linear fitting function. Our fitting function 
applied to this sequence gives acceptable fits. We will try to estimate the error 
in this approach by fitting a curve to data of the form L m ; n < L < 128 with 
L m in = 6, 8, 12, ie for three sets of data. Leaving out more points gives the fitted 
curve an unconvincing look. This gave 

K c = 0.550425 ±0.000025 

which agrees with the previous interval. This estimate is a little lower than that 
of [JV97] (who also provide a nice table of previous results) but their data is 
based on rather small graphs, L < 36. On the other hand, our estimate ends 
up right in the middle of the (rather wide) interval given by GE94J . 

From the interval above we choose the mid-point as our limit, ie we set 
K c = Co = 0.550425, and fit all points to determine the remaining parameters. 
Using the same limit we fitted curves to the K~ (discarding L — 128) and K + 
data. We received the curves shown with the points in the left plot of Figure ITTI 

A different behaviour is expected from the three sequences of energies; U + , 
XJ- and U*. Here U+ -> U+ and U~ -> U~ and the difference [/ ± = U+ - U~ 
should converge to the latent heat U^ 1 , whereas U* should converge to some 
value U c between U~ and U + . Again we see a possibly too big jump in the 
data for U~~ at L = 128 so we will discard this point. Applying the process we 
described above we received 

U~ = 0.5322 ±0.0013 
U+ = 0.5670 ±0.0004 
U c = 0.5513 ±0.0008 
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Again, using the mid-points as limits we received the curves shown with the 
data points in the right plot of Figure [TTJ Using these estimates of 17+ and U~ 
also provides us with estimates of the asymptotic latent heat J7 ± . This resulted 
in 

Ut = U+ - U~ = 0.03480 ± 0.0017 

which is clearly smaller than the estimate 0.0538 (after division by 3) of |J V97j . 
but, as we recall, their estimates were based on much smaller systems. 

0.556 r 




64 32 16 12 8 6 64 32 16 12 8 6 



Figure 11: Left: couplings K , K* , K + (downwards) vs 1/L. Right: energies 
U + , U* , U~ (downwards at y-axis) vs 1/L. 

Applying this procedure to the free energy, where T(K*) — > T Cl and the 
entropy, where S{K*) — > S c , we obtained 

T c = 1.77018 ±0.00005 
S c = 0.86020 ±0.0015 

Note the considerably larger error in S c but also that the value is consistent 
with taking 

S c = T c -ZK c Uc = 0.85983 ± 0.0025 

though we receive a slightly larger error estimate. The data points and the 
fitted curves are shown in Figure fT2l 




64 32 16 12 8 6 64 32 16 12 8 6 



Figure 12: Free energy T(K*) (left) and entropy S(K*) (right) vs 1/L together 
with fitted curves. 

For the magnetisations fi + and fi* we should see a behaviour analogous 
to that of the energies and we have treated them as such. We assume that 
/i~ — > fi~ = 1/3, (i* — » \i c and /i + — » /i+. Continuing with our trusted approach 
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we received 



/Lt+ = 0.3986 ±0.0015 
He = 0.3711 ± 0.0027 

and the fitted curves are shown with the data points in Figure 16.31 
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Figure 13: Left: magnetisations [i* , /i (downwards at y-axis) vs 1/L. 
Right: coefficients for our fitted curves 
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